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Department of Aeronautics and Astronautics, Room 33-311 
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Introduction 


This report summarizes the work done in the period June 1993 through December 1993 on 
the “Modeling of Outgassing and Matrix Decomposition in Carbon-Phenolic Composites” 
program. 

, \ ' 

Progress -/ 

' ' \ 

jjraduate, student David Tai has completed his Master’s Thesis. A new release rate 
equation to model the phase change of water to steam in composite materials was derived 
from the theory of molecular diffusion and equilibrium moisture concentration. The new 
model is dependent on internal pressure, the microstructure of the voids and channels in the 
composite materials, and the diffusion properties of the matrix material. Hence, it is more 
fundamental and accurate than the empirical Arrhenius rate equation currently in use. The 
model was mathematically formalized, and integrated into the thermostructural analysis 
code CHAR. Parametric studies on variation of several parameters have been done. 
Comparisons to Arrhenius and straight-line models show that the new model produces 
physically realistic results under all conditions. 

This work was presented and published at the 1993 JANNAF Rocket Nozzle Technology 
Subcommittee Meeting in October 1993, 1 and was published as a Master’s Thesis a month 
later. 2 The paper and thesis have already attracted wide attention in the rocket nozzle 
industry. Approximately 20‘reprints of the paper, and several copies of the thesis, have 
been requested and distributed. A copy of the paper is attached. The thesis is available 
from the M.I.T. archives, or directly through the TELAC laboratory. 


Current Status 

David Tai has completed his Master's thesis, and has left M.I.T. Graduate student David 
Shia will be taking his place. He has initiated further work in this area, beginning with a 
study of the consequences of non-equilibrium moisture content in the virgin composite. 
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A COMBINED FICKIAN DIFFUSION-DARCY FLOW MODEL 
FOR THE MOVEMENT OF GASES WITHIN A POROUS MEDIA 

David S. Tai, Research Assistant and Hugh L. McManus, Boeing Assistant Professor 
Massachusetts Institute of Technology, Cambridge, MA 02139 

ABSTRACT 

The evaporation of absorbed water has been shown to be a critical factor in the failure of composite 
material ablators. A new reaction rate equation to model the phase change of water to steam in composite 
materials is derived from the theories of molecular diffusion and equilibrium moisture concentration. The 
new model is dependent on internal pressure, the microstructure of the voids and channels in the composite 
materials and the diffusion properties of the material. Hence, it is more fundamental and accurate than the 
empirical equations currently in use. The model and its implementation into the thermostructural analysis 
code CHAR are described. Results of parametric studies on the variation of several parameters are 
presented. 


INTRODUCTION 

When composite nozzle liner materials (geometry shown in Figure 1) are exposed to the high 
temperature environment inside a rocket nozzle, different reaction zones develop as shown in Figure 2. 

First, the material on the heated side will begin to decompose and form a pyrolysis zone. When the 
decomposition finishes, it leaves a layer of char behind. As more heat conducts into the material, the 
pyrolysis zone advances deeper into the virgin material. Ahead of the pyrolysis zone, trapped moisture in 
the virgin material is released. A moisture evaporation zone will also develop and advance ahead of the 
pyrolysis zone in lower temperature material. Gases, which are produced by pyrolysis decomposition and 
moisture evaporation, flow to the surface and provide some cooling, but can cause large internal pressures. 

The thickness of the composite insulation is designed so that the char layer will not reach the back side 
of the material before the rocket engine is shut down. However, several anomalous events can occur during 
the flight which can cause the insulation to fail prematurely. One of the severe anomalies is known as ply- 
lift. Ply-lift refers to the across ply failure of the matrix material and it has been observed in the exit cone 
liners of post-fired rocket engines. The ply-lift failure mode usually occurs in composites with low ply 
angles in the region just underneath the pyrolysis zone. 1 

The ply-lift failure has been attributed to the following mechanisms. When the material is heated 
rapidly, gases are generated and trapped. These excess gases cause a large increase in internal pressure 
which forces the plies apart Since ply-lift usually occurs in carbon-phenolic composite materials at 
temperature below 400 °C and pyrolysis reactions usually do not begin below 400 ®C, it is suspected that 
the high pressure is mainly caused by steam released by absorbed water. When the matrix material in the 
composite decomposes to char, the material's permeability can increase by as much as seven orders of 
magnitude. Thus, the pyrolysis gases inside the pyrolysis zone can escape easily, while steam released in 
the evaporation zone has more difficulty escaping since it has to pass through the relatively impermeable 
virgin material between the pyrolysis and moisture evaporation zones. Large internal pressures are built up 
by gases trapped between these zones, and it is this narrow region where ply-lift failure usually occurs. 
Since the ply-lift failure is caused by high internal pressure, better modeling of the moisture evaporation 
process will result in more accurate prediction of the internal pressure and ply-lift failure. 

BACKGROUND 

In general, a chemical reaction rate can be modeled by an n-th order Arritenius rate equation. The 
Arrhenius rate equation is not dependent on pressure, but the boiling point of water is. To model a 
temperature and pressure dependent moisture evaporation rate equation, McManus 2 proposed using a 
simple straight-line model which reaction rate is constant Another method he proposed is an Arrhenius 
equation with E a as a function of pressure. 3 Arrhenius rate reactions coupled with models for condensing 
of pressurized water have also been proposed. However, all these methods are empirical. There is no 
guarantee that they give an accurate moisture evaporation rate and they provide no insight into rate 
determining mechanisms. So in this study, a new moisture evaporation model based on a more 
fundamental and physical modeling of moisture diffusion and equilibrium moisture concentration is derived 
to improve the accuracy of predicting the moisture evaporation rate. This model will be coupled into an 
existing thermo-chemical-structural analysis program to provide a new and more accurate tool for 
predicting the behavior of ablative composite materials. 
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Then 


R" 1 R' 10 , 2 , „ 

R r R d 9 K 

where R' indicates the derivative of R with respect to r and 9 indicates the time derivative of 9. For the 
R(r) part, 

r 2 R" + rR' + X 2 r 2 R = 0 (( 

This ordinary differential equation (ODE) has the form of Bessel differential function of order zero, 5 and 
the general solution is 

R(r) = Af 0 (lr) + BY 0 (lr) f 

where A and B are constants, J 0 is a Bessel function of order 0 of the first kind and Y 0 is a Bessel function 
of order 0 of the second kind The boundary condition at r = rp requires R(r P ) = 0, and the boundary 
condition at r = r a requires R\r a ) = 0 . These boundary conditions can only be satisfied if the following 
characteristic equation is satisfied (details in Reference 6) 

Jo(Kr p )Y x a H r a )- Y 0 (X n r p )M k n r a ) = 0 (! 

To simplify this characteristic equation, let 


where 


K r P = 


a = r °/ r P 


K r a = az. 


/ 0 (zJF 1 (az n )-r 0 ( Z)t )y l (a z J = 0 (11) 

The roots z* are solved from Eq. 1 1 numerically for any given a , and A, is calculated by Eq. 9. Then Eq. 7 
becomes 

ic^ir) ( 12 ) 

n=l 

where 

y*(r) = Y o a n r p )J Q {X K r) - J 0 (X H r p )Y 0 (X H r) (13) 

R(r) is represented by the sum of the magnitude C„ times the mode shapes y K (r)_ 

For the 9(t) part, 

9 = -XldG (14) 


Substituting Eq. 2 into Eq. 14 


e = -x\d Q exp(--*-)e 

K ! 2 

which is an ODE with solution 


0(f) = 0(0) exp(-doA 2 B g(0) 


where 


* <0 = £ ex *'VM ,,b <17) 

Since the initial conditions is c(0) = c<), and R(r) is assumed to be 1 for r p < r < r a at time zero, 9(0) = 1 and 

0(O = exp(-doA 2 n g(O) (18) 


one term will give good results. We use five terms of Eq. 27, which should give excellent results under 
realistic conditions. 


SURFACE EQUILIBRIUM MOISTURE CONCENTRATION 


In general, the boundary condition at r ■ = r p will not be zero. The equilibrium moisture concentration is 

n hv the following emnirieal eniiaiionA ' 


given by the following empirical equation,' 

( P > Y 

C *" Cm “UJ7’), 


(30) 


where P , is the vapor pressure of water, P sat is the saturated pressure of water, and is the maximum 

moisture content. If we assume that the moisture evaporation rate at the surface is very fast, then 
equilibrium is achieved at once at the surface. We have a new boundary condition 

c(>» = c. (31) 

We set b equal to 1 and equal to c$. Even though supersaturated steam may exists inside the pore 
channel (i.e. P, > P ja /T ) ), c„ cannot be greater than since the material could not physically absorb 
more moisture than the maximum amount In that case, the supersaturated steam will probably condense to 
water inside the pore channels. This possible phenomena is neglected here so that a single phase flow of 
gas can be used. 


MOISTURE EVAPORATION RATE 


The solution for moisture diffusion and concentration derived in the previous section assumed a 
boundary condition of zero concentration at the surface. The solution for realistic boundary conditions can 
be found by a convolution integral (see Reference 6 for details). The degree of conversion is given by 


C„(0 = /(0-J' 

= fit) - £/(«)£ A. exp(-A 2 .doU(0-g(«)])d« 


(32) 


where equilibrium degree of conversion jTO = c. / c 0 is given by 


/ ( , ) = £mii.(A.)» 

P ¥ < Psat 

Psat 

fit) = 

P,2P,at 



(33) 


and /f t) is always constrained to be less than or equal to one. The beta of moisture or degree of dry-out is 
A. = l-C„ (34) 

And the desired stream generation rate is 

dC w 


G w = -c 0 p. 


d? 

dt ' = CoP, ~di 


(35) 


NUMERI CAL M ET HO D 

The steam generation rate (Eq. 35) is calculated numerically by a routine that is embedded in the 
CHAR code. Temperature and pressure conditions are provided by CHAR at each time step. These 
determine the boundary conditions. The degree of conversion C w is found by numerically integrating Eq. 
32 using the same time steps as the rest of the CHAR solution, and the change in degree of conversion from 
the previous time step provides the generation rate. Details of the numerical method are provide in 
Reference 6. 
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moisture release and pyrolysis model. 3 Only the first two parts, which are assumed to involve moisture 
release, were used. The pyrolysis model used was the same in all cases. 

Figure 10 shows the maximum pressure differences calculated for the new, straight-line and Arrhenius 
models. The Arrhenius model overpredicts the maximum pressure difference because it is not pressure 
dependent The straight-line model overpredicts the maximum pressure difference early in the simulation 
because it is not time dependent A straight-line model, especially with AT = 100 K, gives good agreement 
at low temperature rise rates typical in the later part of the simulation. At high heating rates prevalent in the 
early part of the simulation, the straight-line model predicts faster moisture release, resulting in a higher 
maximum pressure difference. 


CONCLUSIONS 

A new method for calculating moisture release rates is based on a micro-scale model of moisture 
diffusion to a nearby pore channel and moisture evaporation on the pore channel surface. The diffusion of 
moisture causes the moisture release rate to be both time and temperature dependent The equilibrium 
condition on the pore channel surface causes the moisture release rate to be dependent on pressure. The 
method has been expressed mathematically and implemented as a module of the CHAR code. It was found 
that only a few terms of the series solution were necessary for accurate results, resulting in good 
computational efficiency. 

The inclusion of diffusion in the model slows the release of moisture. The delay of moisture release 
causes moisture to be released at higher temperatures. This tends to remove the separation between the 
moisture evaporation and char zones, and results in lower predicted pressures. 

The geometry of the pores strongly affects the moisture release rate. Pore spacing has a larger effect 
than the pore size. Larger spacing (or smaller pore size) slows diffusion to the pores and results in lower 
predicted pressures. Very large pore spacing slows the diffusion so much that the effect of moisture on 
predicted pressure is almost lost. Very small pore spacing allows very rapid moisture release to the limit 
that the diffusion effect is lost. 

The value of the diffiisivity rate constants have a lesser effect on the moisture release rate. Varying the 
values of the activation energy E w well outside the measured range had only a moderate effect on the 
moisture release and pressure. Varying the rate constant do by two orders of magnitudes changed the 
calculated pressure difference by a factor of two. 

Comparison of the new model with existing Arrhenius and straight-line models illustrates how the 
more fundamental nature of the new model produces more realistic results. The new model is time, 
temperature, and pressure dependent, and so it produces physically realistic results under all conditions. 
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